#include <kuai/misc/typedef.hpp>

#ifndef _KUAI_SIM_E_FUNCTOR_HPP_
#define _KUAI_SIM_E_FUNCTOR_HPP_


namespace kuai {

	class FunctorHarmonic
	{
	public:
		inline RealNumber operator()(RealNumber r, const RealNumber k[]) const {
			r -= k[0];
			return k[1]*r*r;
		}
		inline RealNumber operator()(RealNumber r, const RealNumber k[], RealNumber& dedr) const {
			r -= k[0];
			dedr = 2 * k[1] * r;
			return k[1] * r * r;
		}
		inline RealNumber operator()(RealNumber r, const RealNumber k[], RealNumber& dedr, RealNumber& d2edr2) const {
			r -= k[0];
			d2edr2 = 2 * k[1];
			dedr = d2edr2 * r;
			return k[1] * r * r;
		}
		enum {
			PAR_NUMBER = 2
		};
	};

	class FunctorQuardric
	{
	public:
		inline RealNumber operator()(RealNumber r, const RealNumber k[]) const {
			r -= k[0];
			return ((k[3]*r + k[2])*r + k[1])* r * r;
		}
		inline RealNumber operator()(RealNumber r, const RealNumber k[], RealNumber& dedr) const {
			r -= k[0];
			dedr = ((4 * k[3] * r + 3 * k[2]) * r + 2 * k[1]) * r;
			return ((k[3]*r + k[2])*r + k[1]) * r * r;
		}
		inline RealNumber operator()(RealNumber r, const RealNumber k[], RealNumber& dedr, RealNumber& d2edr2) const {
			r -= k[0];
			d2edr2 = (12 * k[3] * r + 6 * k[2]) * r + 2 * k[1];
			dedr = ((4 * k[3] * r + 3 * k[2]) * r + 2 * k[1]) * r;
			return ((k[3]*r + k[2])*r + k[1]) * r * r;
		}
		enum {
			PAR_NUMBER = 4
		};
	};

	class FunctorLJ12_6
	{
	public:
		inline RealNumber operator()(RealNumber r, const RealNumber k[]) const {
			RealNumber rr = k[0] / r;
			rr = rr * rr * rr;	// r^(-3)
			rr = rr * rr;		// r^(-6)
			return k[1] * rr * (rr - 2);
		}
		inline RealNumber operator()(RealNumber r, const RealNumber k[], RealNumber& dedr) const {
			RealNumber rr  = 1/r;
			RealNumber rr6 = rr * k[0];
			rr6 = rr6 * rr6 * rr6;	// r^(-3)
			rr6 = rr6 * rr6;		// r^(-6)
			RealNumber krr6 = k[1] * rr6;
			dedr = 12 * krr6 * (1 - rr6) * rr;
			return krr6 * (rr6 - 2);
		}
		inline RealNumber operator()(RealNumber r, const RealNumber k[], RealNumber& dedr, RealNumber& d2edr2) const {
			RealNumber rr  = 1/r;
			RealNumber rr6 = rr * k[0];
			rr6 = rr6 * rr6 * rr6;	// r^(-3)
			rr6 = rr6 * rr6;		// r^(-6)
			dedr = 12 * k[1] * rr6 * (1 - rr6) * rr;
			d2edr2 = 12 * k[1] * rr6 * (13 * rr6 - 7) * rr * rr;
			return k[1] * rr6 * (rr6 - 2);
		}

		enum {
			PAR_NUMBER = 2
		};
	};

	class FunctorLJ_9_6
	{
	public:
		inline RealNumber operator()(RealNumber r, const RealNumber k[]) const {
			RealNumber rr  = 1/r;
			rr *= k[0];
			rr = rr * rr * rr;	// r^(-3)
			return k[1] * rr * rr * (2 * rr - 3);
		}
		inline RealNumber operator()(RealNumber r, const RealNumber k[], RealNumber& dedr) const {
			RealNumber rr  = 1/r;
			RealNumber rr3 = rr * k[0];
			rr3 = rr3 * rr3 * rr3;	// r^(-3)
			RealNumber krr6 = k[1] * rr3 * rr3;
			dedr = 18 * krr6  * (1 - rr3) * rr;
			return krr6  * (2 * rr3 - 3);
		}
		inline RealNumber operator()(RealNumber r, const RealNumber k[], RealNumber& dedr, RealNumber& d2edr2) const {
			RealNumber rr  = 1/r;
			RealNumber rr3 = rr * k[0];
			rr3 = rr3 * rr3 * rr3;	// r^(-3)
			RealNumber krr6 = k[1] * rr3 * rr3;
			d2edr2 = 18 * krr6 * rr * rr * (10*rr3-7);
			dedr = 18 * krr6  * rr * (1 - rr3) ;
			return krr6  * (2 * rr3 - 3);
		}
		enum {
			PAR_NUMBER = 2
		};
	};

	class FunctorPoly5
	{
	public:
		inline RealNumber operator()(RealNumber v, const RealNumber k[]) const {
			return ((((-k[5]*v+k[4])*v-k[3])*v+k[2])*v-k[1])*v+k[0];
		}
		inline RealNumber operator()(RealNumber v, const RealNumber k[], RealNumber& dedr) const {
			dedr = (((-5*k[5]*v+4*k[4])*v-3*k[3])*v+2*k[2])*v-k[1];
			return ((((-k[5]*v+k[4])*v-k[3])*v+k[2])*v-k[1])*v+k[0];
		}
		inline RealNumber operator()(RealNumber v, const RealNumber k[], RealNumber& dedr, RealNumber& d2edr2) const {
			d2edr2 = ((-20*k[5]*v+12*k[4])*v-6*k[3])*v+2*k[2];
			dedr = (((-5*k[5]*v+4*k[4])*v-3*k[3])*v+2*k[2])*v-k[1];
			return ((((-k[5]*v+k[4])*v-k[3])*v+k[2])*v-k[1])*v+k[0];
		}
		enum {
			PAR_NUMBER = 6
		};
	};
}

#endif
